Correlated Levy noise in linear dynamical systems 
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q ■ Linear dynamical systems, driven by a non-white noise which has the 

O \ Levy distribution, are analysed. Noise is modelled by a specific stochastic 

process which is defined by the Langevin equation with a linear force and 
the Levy distributed symmetric white noise. Correlation properties of the 
process are discussed. The Fokker-Planck equation driven by that noise 
is solved. Distributions have the Levy shape and their width, for a given 
time, is smaller than for processes in the white noise limit. Applicability 
of the adiabatic approximation in the case of the linear force is discussed. 
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1. Introduction 

Stochastic dynamical equation (the Langevin equation) describes mo- 
tion of a particle which is subjected to both conservative and stochastic 
force. The latter one can be understood either as a result of elimination 
of internal degrees of freedom or as some external physical process. The 

£-J \ external noise possesses its own time scale and relaxation properties. If re- 

laxation time of processes in the environment is relatively short, the white 

— «: , noise may be a good approximation: the noise variables change rapidly, com- 

pared to the particle variables. Otherwise the Langevin description must 
involve the correlated ('coloured') noise. This problem was widely discussed 
for the Gaussially distributed noise. Well-known physical examples involve 
a phenomenon of narrowing of magnetic resonance lines due to the ther- 
mal fluctuations [1] and the fluctuations of dye laser light [2] . The problem 
of correlated noise also emerges when one eliminates some variables in a 
multi-dimensional dynamical system; then the effective low-dimensional de- 
scription involves correlations even if the original many-dimensional system 
is Markovian [3]. The Langevin equation with the correlated Gaussian noise, 
both additive and multiplicative, is non-Markovian and it resolves itself to 
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an integro-differential Fokker-Planck equation which can be solved exactly 
for simple potentials; otherwise approximate methods may be applied [Sid]. 

Recently, the Levy processes - which constitute a general class of the 
stable processes with the Gaussian process as a special case - attract a 
considerable interest. They are characterised by long tails, which make the 
variance divergent, and can be observed in many systems from various fields: 
porous and disordered materials, hydrology, biology, sociology and finance. 
Realistic problems are usually characterised by high complexity and they 
exhibit collective phenomena; they involve long-range correlations, non-local 
interactions and a complicated, nonhomogeneous (in particular fractal or 
multifractal) structure of the medium. As a result, long jumps may appear 
and the standard central limit theorem is no longer valid. 

It is natural to expect that processes which are driven by a noise with 
long jumps are correlated. As an example can serve an experimental study 
on spontaneous electrical activity of neuronal networks with different sizes 
[5]. It was found that all networks exhibited scale- invariant Levy distri- 
butions. The authors conclude that different-size networks self-organise to 
adjust their activities over many time scales. The power spectrum, calcu- 
lated from the experimental time series, indicates correlations: it obeys a 
power-law decay at low frequencies for all network sizes. 

The non-Markovian master equation governs probability distributions 
in the framework of the continuous time random walk theory [BJ. If jumps 
are Levy distributed, the Fokker-Planck equation is fractional both in time 
and position. The integral operators introduce a competition between sub- 
diffusion and accelerated diffusion; the latter one results from the infinite 
variance. Integral Fokker-Planck equations were solved for both fast and 
slowly decaying memory kernels [7]. They can be generalised to the frac- 
tional orders and to the case of a variable diffusion coefficient [8]. 

In this paper we consider a linear dynamical system which is defined 
by the Langevin equation with the Levy distributed non-white noise. That 
problem was solved by Hanggi and Jung (|3| and references therein) for an 
arbitrary autocorrelation function in the case of the Gaussian noise. How- 
ever, that approach a priori assumes the autocorrelation function and that 
does not exist if a < 2; we will discuss that difficulty in Sec. II. Therefore 
we introduce a specific model of the correlated noise; we require that the 
model process should have the Levy distribution and be correlated (in a 
sense which will be explained in Sec. II). Moreover, it should be as simple as 
possible. We define that process in Sec. II by an adjoint Langevin equation 
which corresponds to the Ornstein-Uhlenbeck process with the white sym- 
metric Levy noise. We also discuss its correlation properties. The Langevin 
equation, driven by that process, is analysed in Sec. Ill for simple forms of 
the potential: the free Levy motion, the constant force and the linear force. 



Results are summarised in Sec.IV. 



2. Ornstein-Uhlenbeck process with Levy noise 

Motion of a particle, which is subjected to the quadratic potential and 
the Levy noise, is described by the following linear Langevin equation 

e(t) = -7£(*) +£(*), (i) 

where the uncorrelated and symmetric noise L(t) is the a— stable Levy pro- 
cess and 7 = const > 0. Eq.([T]), with the initial condition £(0) = 0, can be 
formally solved, 

m= f ' K{t-r)L{dr), (2) 



where K (t) = exp(-jt). The well-known theory of the Brownian motion 
corresponds to the case a = 2. Generalisation to the non-Gaussian stable 
cases, which are defined by Eq.([T]), constitutes the Ornstein-Uhlenbeck- Levy 
process (OULP). If a = 2, trajectories are continuous and Eq.([TJ corre- 
sponds to the standard Fokker-Planck equation. Otherwise jumps emerge 
and their presence requires introducing integral operators. The Fokker- 
Planck equation, which is suited for problems with jumps, contains the 
fractional operator: 

§- t pm=^mm+D-^ P (c,t), ( 3 ) 

where < a < 2 denotes the order parameter of the Levy distribution and 
D > is a constant noise intensity. The Levy distribution itself is given by 
the following Fourier transform: 

1 r°° 
P{L) = - exp(-Dk a )cos{kL)dk. (4) 
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The density distribution p(£,£) can be evaluated either directly from Eq.([2]) 
[9] or by solving Eq.© [10]. The characteristic function reads 



l r°° 
p(fc,t) = — / p(£,i)e- ife ^ = ex P 
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Expression ([5]) corresponds to the Levy stable and symmetric process and 
the width converges with time to a constant, producing a stationary distri- 
bution. The second moment is divergent, unless a = 2, and also the mean 
is divergent if a < 1. 



The Langevin equation driven by the white non-Gaussian noise was 
studied by several authors, both for linear and nonlinear systems |10t [TT1 
[T2J I13j . It was generalised to the asymmetric Levy noise [9] and to the 
multiplicative noise |14l [T5] . OULP was also discussed in Ref . [16] where 
several fractional generalisations were presented. 

Dynamical relation (pQ) introduces a dependence among process values 
£ at different times: the process £(£) possesses memory. For the Gaussian 
case, the autocorrelation function serves as a measure of the memory loss. 
It is defined [I] as the average along a stochastic trajectory: 

G(t)= lim 1 [ T mC(t + r)dt. (6) 
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G{t) can be evaluated as the inverse Fourier transform from the spectral 
function 

S(w)= lim ^ICMI 2 , (7) 

where £(w) stands for the Fourier transform from £(£), by means of the 
Wiener-Khinchin theorem, 

G(r)=T- 1 [S(co)]. (8) 

For the ordinary Ornstein-Uhlenbeck process, which is given by Eq.([T]) with 
a = 2, the stationary autocovariance function G{t) follows directly from 
Eq.©. It assumes the exponential form [3], 

G(t) = -e" 7|r| , (9) 

7 

which corresponds to the Lorentzian shape of S{uS). The correlation time 
t c = I/7 measures the decay rate of G(t). 

Applying the above formalism to the case a < 2 is problematic since 
the variance a 2 = G(0) becomes infinite. To overcome that difficulty, some 
modifications of the standard covariance definition were introduced. One 
can define [T71 [18] the 'codifference' tx,y = cr^ + cry — ffx-Y> where X, Y 
are stable and symmetric processes. For independent X and Y, tx,y = 0; 
codifference resolves itself to the standard covariance if a = 2. On the 
other hand, one can utilise the Poissonian structure of the Levy process 
to introduce an infinite cascade of Poissonian correlation functions which 
correspond to the autocorrelation function [191 ED]- That function depends 
exponentially on time for OULP, Eq.([T]). Standard correlation formalism of 
the general Levy case may be applied if Levy measure in the Levy-Khinchine 
formula [2T] possesses a cut-off |22j : all moments are then finite. Solutions 



of the Langevin equation, which is driven by noise with such a truncated 
distribution, are identical with those for the stable noise up to arbitrarily 
large distances [23] . 

The usual definition of the autocorrelation function, Eq.([6]), may still be 
applicable to the general stable Levy case, despite divergent variance. The 
characteristic function of the increment £(£2) — £(*i) can De formally derived 
[lUj : that function contains all information about two-point correlations. 
Special methods of spectral analysis were developed to handle experimen- 
tal time series which involve long jumps, e.g. calculating the count-based 
periodogram [23]. That method allows one to calculate the autocorrelation 
function and power spectrum for long signals, also containing nonstation- 
ary trends [5j. We will demonstrate, by means of numerical simulation of 
stochastic trajectories, that speed of memory loss for the process (JT]) can be 
determined by means of the ordinary spectral analysis. Let us calculate the 
power spectrum, Eq.flJJ), from a trajectory which follows from Eq.Q and 
has a given length T; the Fourier transform is simultaneously evaluated. 
The relative normalisation of S(u), Sq = 5(0)7 2 , is finite in any calcula- 
tion since T is always finite. However, it depends on T and then cannot be 
determined, as expected. The analysis shows that the quantity S(co)/So is 
well determined in the limit T — > 00, it obeys the Lorentz function 

lim SH/S = l/(7 2 +^ 2 ). (10) 

T->oo 

The renormalised S(oj) is presented in Fig.l for T = 10 and some values 
of a and 7. All curves follow the Lorentzian shape. The value of So, which 
emerges from that calculation, may be large, it ranges from 1 (a = 2) to 
10 3 (a = 1.2). 

Equivalence of the expression (jSJ) with the ensemble averaged covariance 
is not obvious since a system with long jumps may be non-ergodic [25J. The 
latter quantity can be directly evaluated if one introduces a cut-off in the 
distribution Q. We define the ensemble-averaged autocorrelation function 

C(r) = (£(0)£(r))/(e(0) 2 ) (11) 

on the assumption that P(L) = for L > L c . Fig.l presents that quantity; 
it was derived from the time evolution of individual trajectories by averag- 
ing over the ensemble. The figure demonstrates that also C{t) obeys the 
exponential dependence fl§J). 

3. Langevin equation with coloured noise 

In this section we study the stochastic dynamics of a particle which is 
subjected to the Levy correlated noise and the linear deterministic force. 



C/3 




Fig. 1. Renormalised spectral function for OULP, Eq.((T|), calculated from evolution 
of a trajectory up to t = 10 4 , for the following cases: a = 1.2 (dashed line), a = 1.5 
(green dots) and a — 2 (blue dashed-dotted line). Red solid line denotes the Lorentz 
function (1101) . Upper and lower curves correspond to 7 = 1 and 2, respectively. 
Inset: C(t), calculated from an ensemble of 10 6 trajectories with L c = 10 4 , for 
7 = 1 and 2 (solid lines). Red dashed lines represent the function e _7T . 



The noise £(£) is represented by OULP, Eq.([T]). Then we have to solve a set 
of two Langevin equations, 



x{t) = /o - Xx(t) + 7 f(t) 



(12) 



where 7 > 0, A > and /o are constants. In the presence of jumps, the 
system remains far from the thermal equilibrium and the detailed balance is 
violated. Then £(£) can be regarded as an external noise which has its own 
time scale, determined by the parameter 7. In general, processes which obey 
Langevin equation with the correlated noise are non-Markovian since the 
process values are evaluated from mutually dependent noise increments [3]. 
For large 7 (short correlation time), £ is a fast, rapidly relaxing variable and 
the process can be approximated by a corresponding white-noise problem, 
by using the methods of adiabatic elimination of fast variables [3J H] . 



3.1. The case without deterministic force and with a constant force 

The force-free motion, with the white Levy noise, is a generalisation of 
the Wiener process; it describes simple diffusion if a = 2. Generalisation 
to the coloured noise is defined by Eq, (fT2]) with /o = A = 0. We assume 
the initial conditions x(0) = £(0) = 0. Our aim is to find the probability 
distribution of the variable x. One can solve Eq. (|12|) and utilise the fact 
that x(t) is still a process with independent increments, though multiplied by 
some function of time; then convolution of densities can be performed. That 
method was applied in Ref. [llj to the second order Langevin equation for the 
case a = 1. We apply van Kampen's method of compound master equations 
|13j which consists in solving the joint fractional Fokker-Planck equation for 
the two-dimensional system, (x, £), and integrating over the internal noise £. 
That method is relatively simple in the case without potential and formally 
applicable also to nonlinear systems with a multiplicative noise. In the linear 
case, the existence, uniqueness and positiveness of the solution is ensured 



The Langevin equations (|12p correspond to the fractional Fokker-Planck 
equation for a joint probability distribution p(x,£,t) [26, 14|: 
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-l^ri + l^i + D- 



P(x,£,t). (13) 
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Knowing the solution of Eq. (fT3|) , the probability distribution of the variable 
x can be obtained by integration over all possible realisations of the noise 

/■oo 

p(x,t)= / p(x,Z,t)d£. (14) 



Fourier transformation of Eq. (|13p . in respect to both x and £, produces the 
equation for the characteristic function p(k, k, t), 

d „ d ,_ 

— p-j(k-K)— p = -D\K\ a p, (15) 

which can be solved exactly by the method of characteristics; details are 
presented in Appendix. The Fourier transform of the solution, Eq. (|14p . 
follows from Eq. (|A6p : 



p(k, t) = p(k, 0, t) = e - D ^)\ k \ a , (16) 

where 

1 f g K a , 
a(t) = - / dK (17) 

W 7 h 1 - « 
and g = 1 — e -7 *. Eq. (|16|) predicts the Levy shape with the order parameter 
a. The width parameter a(t) can be estimated in the limit jt S> 1, when 




Fig. 2. Probability distributions at t = 1 for the force-free case calculated by the 
Monte-Carlo simulations (points) for 7 = 1,2, 5, 20 (from top to bottom); the most 
diffused case corresponds to the white noise limit (7 = 00). Analytical results, 
calculated from Eg . (|2"0")) with a from Eg. ([TTJ)) . are presented as solid lines. The 
order parameter a = 1.5. Numerical simulations were performed with the time 
step r = 0.005 and averaged over 10 7 events. 



the main contribution to the integral comes from the vicinity of the upper 
integration limit, since then the denominator is close to zero: 



1 



(T(i)« -(1-e 
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-7t\a 



(18) 



In the limit "ft —> 00, a rises linearly with time and p(x, t) coincides with the 
solution of the uncorrelated problem. Convergence to that solution depends 
on a: it is faster for smaller a. 

The integral (|17p can be exactly evaluated if a is a rational number. In 
particular, for a = 3/2 it yields 



a(t) = - 

7 



:i 



-lt\l/2 



(1 - e" 7 *) 3/2 + arctanh\/l - e~^ 



(19) 



In the limit "ft 3> 1, the expression (|19p predicts a time shift, in respect to 
the white noise case, since it can be approximated by a ~ t— (8/3— 2 In 2) /j. 



Numerical values of the probability distribution p(x, t), which result from 
inversion of the characteristic function (|16|) . can be obtained from the series 



expansion 



° r[l + (2n + l)/a] ' - '-" 



n=0 v ' 

if |x| is not too large. Fig. 2 presents those distributions for the case a = 1.5 
at £ = 1, a(t) was calculated from Eq. (|19|) . Figure shows that the memory 
affects the rate of spreading of the distribution: p(x, t) is broadest for the 
white noise case, 7 = 00, and it contracts to the delta function in the 
limit 7 — > 0. Results are compared with the Monte Carlo simulations of 
individual trajectories, according to the stochastic equations ([12]) . For that 
purpose, a simple Euler algorithm was applied. The white noise value at 
i— th integration step, Li, was represented by the term r 1 ' a Lj, where r was 
the step size [28]. Probability distributions were obtained by averaging over 
an statistical ensemble of the individual trajectories. Since the analytical 
result does not contain any approximation, agreement with the simulations 
is exact. 

Problem of the linear potential, —fox, where /o =const., can be reduced 
to the force-free case which was discussed above. The first equation in 
Eq. (|12p takes the form x(t) = /o + 7^(0- From the corresponding fractional 
Fokker-Planck equation, 



9 r <■ n 
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we derive equation for the characteristic function: 



p(x,t,t), (21) 



^- t p-j(.k-K)-?-p = -(if k + D\K\ a )p. (22) 

Its solution, p(k,K,t) = e~ l f° kt pQ, where po is given by Eq. (|A6p . follows 
from the general theory [J3]. It can be also obtained by separation of real 
and imaginary parts of p(k, k, t) and by solving the resulting set of two 
equations. Integration over the variable £ produces the final result: 

p(k,t)=e~ ifokt p , (23) 

where po(k,t) follows from Eq.(|16p. The distribution p(x,t) has the same 
shape, for any time, as that for the case /o = but it is shifted by fat. That 
means that the average rises linearly with time, (£) = fot (if a > 1), and 
the distribution widens with time according to the function cr(t), Eq. (|17p . 
In the limit 7 — > 0, po(x,t) = 5(x) which corresponds to a deterministic 
motion with velocity /o- Probability distributions which follow from the 
Monte Carlo simulations (not presented) agree with the solution 
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Fig. 3. Exemplary stochastic trajectories in the space (£, x), calculated from Ea. (TT^|) 
with time step r = 5 • 10~ 4 up to t = 3, for A = 1 and 7 = 1. The trajectory for 
the case a = 1.5 is positioned in upper-right quarter of the figure. 



In the limit jt — > 00, Eq.(J23J) coincides with the solution of fractional 
Fokker-Planck equation with the constant force for the white noise case [10] . 
The problem of transport in an effective constant force field emerges in the 
framework of the continuous time random walk theory when one considers a 
biased walk [29] . It resolves itself to the fractional Fokker-Planck equation 
with a drift term. 



3. 2. Linear force 

The system is defined by Eq. (fT2j) with /o = 0, where A > measures 
intensity of the deterministic force. The aim of this section is a comparison 
of exact probability distributions, obtained by numerical simulation of two- 
dimensional stochastic trajectories from Eq. (|12|) . with predictions of the 
adiabatic approximation. 

Fig. 3 presents examples of stochastic trajectories for two cases: the Levy 
distribution with a = 1.5 and for the normal distribution. In the former 
case, large jumps, typical for the Levy processes, are visible along the hori- 
zontal direction which represents OULP (Eq.([T])). The process x(t), in turn, 
is stronger localised for both values of a. The plot shrinks in the horizontal 
direction with increasing 7 (not shown) which reflects the fact that £ be- 
comes the fast variable: it relaxes rapidly to £ = 0. Averaging over a large 
number of trajectories produces the probability distribution p(x,t). Fig.4 
demonstrates that it converges with time to the stationary distribution, as 
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Fig. 4. Time evolution of the probability distribution for the system with linear 
force, Eq. (fT2|) . calculated for the following times: 1, 2, 3, 5 (black solid lines from 
top to bottom). The case t = 10, which corresponds to the stationary solution, is 
marked by red solid line. The stationary solution which is predicted by the adia- 
batic approximation, Eg. (1251) . is shown as blue dashed line. The other parameters: 
a = 1.5, A = 1 and 7 = 1. 



in the white noise case. The time which is needed to reach the steady state 
equals 5 for the case presented in the figure. The shape of p(x,t) coincides 
with the Levy distribution for any 7 and its order parameter a corresponds 
to that of the driving noise L{t). The apparent width rises with 7 and, for 
large 7, the white-noise limit is reached. 

To estimate the dependence (7(7) the characteristic function exp(— a{t)\k\ c 
was evaluated. Results are presented in Fig. 5. The distribution very slowly 
converges with 7 to the white-noise value whereas it shrinks to the delta 
function for 7 — > 0. 

The adiabatic approximation in the case of the normally distributed 
noise was discussed in Ref.[30]; we apply a similar procedure. Combination 
of equations (fl~2]) yields a single second order stochastic equation: 



x(t) = -(A + j)x(t) - Xjx(t) + 7 L(t). 



(24) 



One can demonstrate, by introducing a new time variable t' = yfyt, that 
the term x is small both for 7 — > and 00. Therefore, Eq. (|24|) can be 
approximated by the following equation 



x(t) = —Xc-yx(t) + CjL(t), 



(25) 
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Fig. 5. Width parameter a, evaluated from the characteristic function for t = 1, as a 
function of memory parameter 7 (points). Results of the adiabatic approximation, 
Eq. (l2l))) . are marked by stars. The parameters are: a = 1.5 and A = 1. Horizontal 
line marks the white noise limit. 



where c 7 = 1/(1 + A/7). The corresponding fractional Fokker-Planck equa- 
tion is analogous to Eq.([3]) and it can be easily solved. Fourier transform of 
the solution is p a (k,t) = exp(— a a (t)\k\ a ), where the apparent width 



ffaW 



c a D 



-a\t 



)■ 



(26) 



The adiabatic solution, p a (x,t), converges with time to the steady state 
and it coincides with the uncorrelated process in the limit 7 — > 00; Eq. (|26p 
implies that a a rises with 7. Eq. (|25p is exact both for 7 — >• - when 
the delta function is the solution - and in the limit 7 — > 00 (the Smolu- 
chowski limit). For intermediate values of 7, one can expect that Eq, (j25p 
is a good approximation on time scales t > 1/(A + 7) and at distances 

»D- 1 /2/( 7 V2+A 7 - 1 /2) [31 . 

The width parameter a(t) for the exact solution is compared with a a , 
predicted by Eq. ([26j) . in Fig. 5. Some differences are visible but qualitative 
agreement of the functions (7(7) for both cases is good in the entire range 
of presented 7 values. In general, however, discrepancies may be more 
pronounced. For example, the adiabatic approximation underestimates the 
width of the steady-state distribution for 7 = 1, which is shown in Fig. 4, 
by a factor of two (0.24 vs. 0.48). 
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4. Summary and conclusions 

We have studied the linear dynamical systems which are driven by the 
additive, non-white Levy noise. That noise is modelled by a concrete, sim- 
ple stochastic process, OULP. Then the system is defined in terms of two 
Langevin equations. OULP reveals the memory effects, as for the ordinary 
Ornstein-Uhlenbeck process, but their quantitative description is more dif- 
ficult because of the divergent variance. We have presented a numerical ex- 
ample which demonstrates that the renormalised autocorrelation function 
G(t) can be useful as a measure of the memory loss; it falls exponentially 
with time for any order parameter a. The same result was obtained for 
the ensemble-averaged autocorrelation function on the assumption that the 
Levy distribution is truncated. 

In the absence of any deterministic force, the non-Markovian problem 
resolves itself to the Wiener-Levy process (correlated Levy motion). The 
resulting probability distribution has the Levy shape, with parameter a, and 
it converges with time to that for the uncorrelated case. Correlation time 
r c = I/7 determines the distribution width: the larger r c , the narrower the 
distribution. The case of the constant force /o is similar; shape and width 
of the distribution is the same but the time-dependent shift fot emerges. 

Solution for the case of the linear force converges with time to the steady 
state, as for the white-noise problem, and its shape is Levy with parameter 
a. Inclusion the finite correlation time narrows the distribution, analogously 
to the case without a force. The above observations agree with the adiabatic 
approximation approach. That method deals with a corresponding, effective 
white-noise process and resolves itself to the Langevin equation of the first 
order. It is supposed to be accurate if 7 is sufficiently large or if 7 — ^ 0. For 
intermediate values of 7, overall predictions of the adiabatic approximation 
in respect to the distribution shape and its dependence on 7 are still correct, 
nevertheless some quantitative discrepancies have been found. 

APPENDIX 

In the Appendix, we solve the fractional Fokker-Planck equation, Eq. (|15|) . 
by means of the method of characteristics. 
First, we put the equation into the form 

\K\~ a —p(k, K, t) - j(k - K)\K\- a —p(k, K, t) = -Dp(k, K, t). (Al) 



Eq. (jAip is the linear partial differential equation of the first order with only 
two variables, t and k, since k can be regarded as a constant parameter. 
The equation can be handled by the method of characteristics [3T]. The 
method consists in reducing the problem to solution of a system of ordinary 
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differential equations (characteristic equations) . Those equations determine 
variables t, £ and z, as functions of parameters s and r, on a characteristic 
curve. They are of the form 

dt(r, s) __ 

ds 
dn(r, s) 

ds 

dz(r, s) 

ds 

with the initial conditions 



-~f{k - k)\k\-° (A2) 

-Dz 



t(r,0) = 

«(r,0) = r (A3) 

z(r,0) = 1; 

the third condition reflects the requirement that p(x,£, 0) is to be the delta 
function in the variable £. We must solve the system ()A2j) and then eliminate 
the parameters r(£, k) and s(t, k). The final solution of Eq. (|Alj) is given by 
p(k,K,t) = z(r,s). Combination of the first and second equation gives the 
relation between t and k on the characteristic curve: t = ln[(/« — k)/(r — 
k)]/j, where the initial conditions (|A3[> were taken into account. The above 
relation determines the parameter r: 

r(t,K) = k-(k- K )e~~ ft . (A4) 

Integration of the third equation (|A2p is straightforward, z(r, s) = e~ Ds , and 
s, as a function of the variables k and t, follows from the second equation: 

s(t,K) = - ' ' , dK ; . (A5) 

7 J r k' - fe 

The final solution reads 

p(k, K ,t) =e~ Ds , (A6) 



where s is given by Eq. (|A5p . The solution (|A6|) can be verified by a direct 
inserting into Eq. ljAip and applying the Leibniz rule for differentiation of 
the integral. 
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